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Multiple avalanches, initiated by simultaneously toppling neighbouring sites, are studied in three 
different directed sandpile models. It is argued that, while the single avalanche exponents are 
different for the three models, a suitably defined two-avalanche distribution has identical exponents. 
The origin of this universality is traced to particle conservation. 

PACS numbers: 05.65.+b,05.70.Ln,45.70.Ht,47.27.eb 



The sandpile model is a paradigm for self-organized 
criticality wherein long range correlations are generated 
without any parameter being fine tuned P, Q . The orig- 
inal version of the model and its variants (see for a re- 
view) have a common feature: slow driving during which 
particles are added to the system, and fast dissipation 
during which the system relaxes through avalanches. The 
steady state is characterized by power law correlations. 

Conservation laws are known to constrain correlation 
functions of driven-dissipative systems. A well known 
example is the Kolmogorov 4/5-th law of three dimen- 
sional fluid turbulence 0, S S 01 • The conserved quan- 
tity is energy which is pumped in at large length scales 
and dissipated through viscosity at small length scales. 
The 4/5-th law states that, in the inertial range [dis- 
tances r between driving and dissipation length scales], 
{[vi{r,t) — vi{{),tyf') = — fe'^, where vi(f,t) is the longitu- 
dinal component of the velocity at point r at time t, and 
e is the energy dissipation rate. The linear dependence 
on r remains true in all dimensions, while the proportion- 
ality constant is a function of dimension. There are other 
examples, mainly from turbulence, of a correlation func- 
tion being determined by the constant flux of a conserved 
quantity. Examples include magneto- hydrodynamics 
burgers turbulence Q and advection of a passive scalar 
(see and references within). These relations are cen- 
tral to understanding turbulence, acting as checkpoints 
for phenomological theories. Examples outside turbu- 
lence are few. In a recent paper [10|, this relation was 
generalised to an arbitrary driven dissipative system that 
showed the general features of turbulence. Exact results 
were obtained for specific models, namely wave turbu- 
lenc e and models of diffusing-aggregating particles 



In sandpile models, the total number of particles is con- 
served in each toppling. As a consequence, can any corre- 
lation function be determined? In this paper, we answer 
this question in the context of directed sandpile models. 
Consider multiple avalanches obtained by adding parti- 
cles simultaneously at nearby lattice sites. It is argued 
that a suitably defined two avalanche joint probability 
distribution function [defined later] will play the role of 
the three point velocity correlations in the Kolmogorov 
4/5-th law, and will have a scaling exponent which is in- 



dependent of dimension and hence identical to the mean 
field answer. 

We define the three sandpile models studied in this 
paper on a directed square lattice of horizontal extent 
L and vertical extent T (see Fig. [Ij. Periodic bound- 
ary conditions are imposed in the x-direction and open 
boundary conditions along the ^-direction, also referred 
to as the time direction. The number of particles at a 
site {x,t) is denoted by a non-negative integer h(x,t). 
All the three models are driven by adding a particle to a 
randomly chosen site on the top layer {t — 0) and then 
letting the system relax according to the following rules 
of evolution. 

The deterministic model 14]: A stable configuration has 
all h{x, t) = 0, 1. If (/i, x, t) > 2, then it relaxes by trans- 
ferring two particles, one each to its two downward neigh- 
bours, i.e, h(x, t) decreases by 2 and h{x ^ l,t + 1) and 
h{x -f 1, t -I- 1) increas e by one. 

The stochastic model 15, 1^: This model has the same 
rules of evolution as the deterministic model except for 
one difference. The toppling is now stochastic. When a 
site (x,t) topples, with probability 1/4 both particles go 
to (x — 1, i -f 1), with probability 1/4 both particles go 
to {x + l,t + 1), with probability 1/2, {x — l,t + 1) and 
(x + 1, t + 1) receive one particle each. 
The sticky model fl7i ]: In this model, the heights can 
take any non-negative integer value. A site is consid- 
ered unstable if h{x,t) > 2 and it received at least one 
particle the previous time step. All unstable sites re- 
lax simultaneously as follows: With probability p, the 
height decreases by 2 and a particle is added to each of 
its downward neighbours. With probability (1 — p), the 
site becomes stable without losing any particles. 

In all the three models, if a site at the bottom {t = T) 
topples, then the height at that site reduces by two, 
and the two particles are removed from the system. An 
avalanche is defined as the number of topplings that 
the system undergoes after a particle is added to a sta- 
ble configuration. In the steady state, the probability 
of an avalanche of size s is a power law distribution 
P{s,T) ~ s~''f{sT-^), when L > Tl/^ where z is 
the dynamic exponent. These two exponents are not in- 
dependent from each other. Particle conservation from 
layer to layer results in the scaling relation (s) ~ T (for 
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FIG. 1: The directed square lattice with L sites along the 
2;-axis and T sites along the t-axis. Periodic boundary condi- 
tions are applied in the horizontal direction. 



example, see Q), implying that (5(2 ~ t) — 1. 

The three models belong to three different universal- 
ity classes. For the deterministic model, first studied in 
Ref. r = 4/3, 5 = 3/2 in d = 2, r = 3/2, S = 2 
in d > 3. In d = 3, the mean field results have logarith- 
mic corrections . Stochasticity in the toppling rules is 
known to change the universality class of sandpile mod- 
els For the stochastic model, it was argued that 
T = 10/7, 5 = 7/4 in d = 2, T = 3/2, 5 = 2 in d > 3, with 
d = 3 having the mean field exponents with logarithmic 
corrections [ll, [3] • The sticky model was introduced in 
Ref. [r?! . Introducing stickiness changes the universality 
class of the sandpile model away from deterministic and 
stochastic classes [3| ■ The avalanche exponents are then 
related to the exponents of directed percolation. From 
the best numerical estimates for directed percolation ex- 
ponents, it was shown that r « 1.32, 6 — 1.47 in d = 2 
[17| . In addition to having different exponents, the sticky 
model is not abelian, unlike the other two models. 

We now define the two-avalanche distributions. Con- 
sider avalanches initiated by adding two particles simul- 
taneously at nearby lattice sites (denoted by 1 and 2) 
on the top level. Let the set of sites belonging to the 
avalanche associated with site 1 (site 2) be denoted by 
Si ('S'2). When a site topples, if it had received particles 
from only sites in Si (S'2), then the site is assigned to 
(S'2). If on the other hand, it had received particles 
from sites belonging to Si as well as S'2, then the site 
is assigned randomly to one of the sets. Let si (S2) de- 
note the number of topplings undergone by sites in Si 
(S'2). We will denote the joint probability distribution 
by Pi(si,S2). For abelian models, we can also define a 
two avalanche distribution as follows. Topple a site. Let 
the avalanche size be si. Then topple the neighbouring 
site. Let the avalanche distribution be S2- Let the joint 
probability distribution be denoted by P2(si, 52)- In this 
paper, it is argued that Pi and P2 have scaling exponents 
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FIG. 2: (a) The avalanche when site A is toppled is shown 
by filled circles (black or grey). The black circles have height 
after the avalanche. The grey circles have the same height 
as before the avalanche. The circles with diameter drawn are 
not part of the avalanche, but will have height 1 after the 
avalanche, (b) The avalanche when B is toppled after A is 
shown by hatched circles. The left boundary of the avalanche 
B is adjacent to the right boundary of avalanche A. 



3 in all dimensions, i.e., 

P^iAsl,As2)^^P^{Sl,S2), 1 = 1,2. (1) 

We give a heuristic argument supporting this conjec- 
ture. Consider the single avalanche probability P{s,t). 
In continuous time, it schematically obeys the equation 



dP{s,t) 
dt 



j dsids2Pl{si,S2)5{si + S2 - s). 



(2) 



where Pi is the two avalanche distribution defined above. 
Use the fact that, for all the three models (s) ~ Jt, where 
J is a constant [3|. Multiply Eq. ([2]) by s and integrate 
over s. The left hand side is a constant independent of t 
and s. Equation ([2]) then reduces to 



const ~ j dsdsids2sPi{si, S2)S{ 



Sl + S2 - s). (3) 



A dimensional analysis of the right hand side of Eq. ([3]) 
immediately predicts 



Pi[s,s) T, 



(4) 



or more generally Eq. ([T]) with z = 1. For the abelian 
versions of the model, the order of toppling is not crucial. 
Hence, one can conjecture that instead of simultaneous 
toppling, the toppling could be sequential and that P2 
obeys the same scaling law as in Eq ([T]). 

We now give a direct proof that, in the deterministic 
model, ^2(51,82) has the scaling as in Eq. ([T]). In the 
steady state, each configuration of the model has equal 
weight [ll]. Thus, each height is and 1 with probabil- 
ity 1/2 independent of other sites. The avalanches then 
have no holes i.e., an avalanche is described by the two 
boundaries, each of which are random walkers that an- 
nihilate on contact. The clusters also have the following 
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FIG. 3: The variation of (a)Pi(s,Si) and (b) P2(s,Si) with s 
is shown for the deterministic model for si — s,s%/2, s- The 
bottom two curves have been shifted for clarity. The solid 
hnes are power laws with exponent (a) 2.97 and (b) 2.98. 




FIG. 4: The variation of (a)Pi(s,si) and (b) P2{s,Si) with 
s is shown for the stochastic model for si — s,s\/2, s- The 
bottom two curves have been shifted for clarity. The solid 
hnes are power laws with exponent (a) 3.03 and (b) 3.00. 



property. Consider the right boundary. If it is at {x,t), 
the next time step, it can either go to {x + l,t + 1) or 
{x-l,t+ 1). If it goes to (x + 1, < + 1), then h{x, t) = 0. 
It it goes to {x— 1, t+1), then h{x + \, t+1) — 1. Similar 
rules exist for the left walker. An example is shown in 
Fig. mja) with site A having been toppled. Now consider 
the case when B (see Fig. [UJb)) is toppled. The black 
circles cannot topple because they have height 0, ensur- 
ing that the two avalanches do not overlap. On the other 
hand, the sites with height 1 will necessarily topple pro- 
vided the avalanche survives up to that level. Hence the 
right boundary of first avalanche and the left boundary 
of second avalanche will be adjacent to each other (see 
Fig.Hb)). 

The calculation of P2(si, S2) now reduces to the prob- 
lem of three annihilating walkers. Let us calculate the 
probability that both avalanches exceed time t. This 
is equal to the survival probability of three annihilating 
random walkers up to time which varies as t^^^'^ when 
t > 1 ;2Q]. Using the scaling s ~ ^3/2 J^J]^ 

we obtain that 

ST ir dsids2P2isi,S2) s~^, or P2{s, s) ^ s^^, consis- 
tent with Eq. (III). The argument for Pi(si,S2) proceeds 
on exactly the same lines and we omit the argument here. 

For the other two models, we rely on Monte Carlo 
simulations. Simulations were done for a lattice with 
L — 1024 and T — 8192. Logarithmic binning was used 
with bin size ln(-\/2). In the steady state, the data was av- 
eraged over 2 x 10* avalanches initiated by toppling near- 
est neighbours. These multi avalanches were interspersed 
with single site avalanches. Avalanches that reached the 
boundary were omitted from the statistics in order to 
prevent strong finite size corrections 21]. 

We study the variation with s of the probability dis- 
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FIG. 5: The variation of Pi(s,si) with s is shown for the 
sticky model for si — s,s\/2, s. The simulations are for p = 
0.85. The bottom two curves have been shifted for clarity. 
The solid lines are power laws with exponent 3.03. 

tributions Pi(s,si) and P2{s,si) with si = s,s\/2,2s. 
The exponents are determined through the maximum 
likehhood estimator method [H, [2^. Let Pi{si,s) ~ 
(sis)""^'/^, for i — 1,2. The numerical estimates for ti 
and T2 are shown in TablelH The data is shown in Figs. [3] 
[deterministic] 0] [stochastic] , and [5] [sticky] , all in good 
agreement with Eq. ([T]). 

In dimensions greater than the upper critical dimen- 
sion, we expect the scaling in Eq. ^ to hold, the mean- 
field avalanche exponent being 3/2. The deviation from 
meanfield should be most pronounced in two dimen- 
sions for which Eq. ([T]) has been numerically verified. In 
other dimensions, we expect that an equation of the form 
Eq. (12) should hold, maybe with a different joint proba- 
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TABLE I: Numerically obtained values for the exponents ri 
and T2 for the different models. 



Model 




T2 


Deterministic 


2.97 ±0.04 


2.98 ± 0.04 


Stochastic 


3.03 ± 0.05 


3.00 ± 0.04 


Sticky 


3.03 ±0.06 





bility distribution. For example, in three dimensions, it 
will be Pi(si, S2, S3). However, the main contribution to 
this three point function will be when one of the s^'s is 
small and we retrieve an effective two-point function. 

To summarize, the two-avalanche distribution was 
studied for three directed sandpile models. While the 
three models have different exponents for the single site 
avalanche distribution, it was shown numerically and 
through a heuristic argument that the two avalanche dis- 
tribution is the same for all three. Exact results were 
obtained for the deterministic model. The robustness of 
the result is due to particle conservation layer by layer, 
leading to the scaling relation (s) ^ and is not depen- 
dent on the details of the model. 

I thank O. Zaboronski, C. Connaughton and D. Dhar 
for helpful discussions. 
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